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We show that in a broad class of directed abelian sandpile 
models that had been expected to have the same exponents 
as the Dhar-Ramaswamy model, the avalanche exponent de- 
pends upon the details of the interaction, calling into question 
the general existence of universality classes in self organized 
critical models. 

64.60.Lx,05.70.Jk 



The existence of universality classes in systems exhibit- 
ing self-organized critical behaviour(SOC) is expected on 
the basis of the similarity of the SOC state to the criti- 
cal state of equihbrium systems. The particular physical 
constraints that might determine these classes are, how- 
ever, very much an open question, but they are thought 
to resemble the universality classes of equilibrium sys- 
tems, at least as far as symmettry and dimensionality 
are concerned. In one dimension, it was shown recently 
by Paczuski and Boettcher ^ that a model for rice piles 
and another for interface depinning were in the same uni- 
versality class, and they conjectured that the Burridge- 
Knopoff model ||^ was contained there as well. Sub- 
sequently, Cule and Hwa |^ described a related model 
that was in the same universality class. A classification 
scheme for 2-dimensional abelian sandpile models was 
suggested by Ben-Hur and Biham In their paper, 
models were classified according to the "sand current" 
,J, which resulted when an unstable site was toppled. 
Models with J = and J were said to be in the 
" nondirected" and " directed" universality classes respec- 
tively, and a distinction was made between those that 
were nondirected only on average. Dhar and Ramaswamy 
l^(DR) had earlier presented an exact solution for a di- 
rected model in two dimensions, and had shown numeri- 
cally that a closely related directed model had the same 
exponents. Ben-hur and Biham described having run a 
third model with the same exponents. However, in their 
simulations, only models with nearest neighbor interac- 
tions were considered. Our investigations of longer (but 
still finite) range interactions in directed models show 
a continuous variation of exponent values ranging from 
those of the exactly solvable nearest neighbor model of 
DR to that of an exactly solvable infinite range interac- 
tion model that we introduce here. 

No simple classification is apparent in these models, 
with exponents depending on the details of the interac- 
tion, as well as the range. While it is not unexpected 
that universality classes in SOC would be more complex 



than those in equilibrium models, this finding calls into 
question the idea of using universality ideas to describe 
these models at all, at least in two dimensions. 

Since the variation in exponents is not large(2.33-2.5), 
we introduce a simulation method which exploits the 
properties of a subclass of the directed models to greatly 
reduce the deviation from power law behaviour seen in 
the simulations as a result of finite system sizes. This 
allows more accurate determination of scaling exponents 
from relatively small system sizes. We emphasize, how- 
ever, that the standard methods for obtaining the expo- 
nents yield the same results. 

To fix notation, we define a sandpile model by three 
variable aspects. 

lattice: In two dimensions the model dynamics occur on 
a MxN lattice CmM composed of d = MN sites 
where sand can build up. This lattice is embedded 
in an infinite lattice of sites which absorb all sand 
which fall on them. The state of the system at 
any given time is a c? dimensional integer valued 
vector whose components represent the height of 
individual sites. 

critical height: The critical height he is the maximum 
amount of sand which can sit on a site and not 
topple. If a site has more than he grains it is said 
to be unstable and will topple. 

toppling rule: The toppling (or relaxation) rule deter- 
mines how sand is removed from an unstable site 
i and redistributed to sites in some predetermined 
neighborhood ,A/i, of i. 

One then defines the set of toppling vectors, {Tj} ,where 
the subscript corresponds to one of the d lattice sites. 
The unstable site i with hc + \ ov more grains is toppled 
by subtracting the vector Ti from the state vector. The 
s*'' component of Ti is defined as 



= ''rr' (1) 

a positive integer if s e A/^ 



where 



his) = 







otherwise 



For two sites i and j, Afj is simply a translated copy 
of Afi. That is, all unstable sites topple identically in the 
sense that their contents are distributed in the same pat- 
tern relative to the toppling sites. The models considered 



1 



here are further constrained so that X]s=i fii^) — 
where the equaUty holds when A/i G CmN'- 
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FIG. 1. Toppling patterns and avalanche exponents for the 
models studied. The exponents are accurate to at least ± .01. 

In our simulations we considered a variety of directed 
models with differing toppling vectors some of which are 
described in Fig.l. Model I was previously solved exactly 
by Dhar and Ramaswamy(DR) 

Abelian sandpile models (ASM) are conventionally 
driven by adding a grain of sand to a random site in a sta- 
ble configuration. If this site becomes unstable then the 
appropriate toppling vector is subtracted from the new 
state perhaps inducing other sites to become unstable. 
On the next time step the lattice is reexamined for un- 
stable sites which are then toppled. The number of time 
steps before no new unstable sites are found is said to 
be the avalanche's lifetime. The total number of toppled 
sites is called the avalanche's size or mass. Once the sys- 
tem relaxes, another grain is randomly added to the new 
stable state. DR showed that once an ASM is driven into 
the stationary state by this method, all reccurent config- 
urations occur with equal probability. So in theory at 
least, one could get the same statistics by adding sand to 
random members of the recurrent set. Typically this is 
no small task since for the general ASM the recurrent set 
is remarkably complex. However, there is an infinite set 
of relaxation rules falling under the directed model clas- 
sification which posess trivial recurrent sets. It can be 
shown that if there exists an orthogonal transformation 
which leaves the toppling matrix, Ajj = [Tj]i, triangular, 
the recurrent set will consist of all possible stable config- 
urations. [0J|]. All of the models we simulated posess this 
trivial stationary state. 

Since an avalanche may be seeded at a site near a 
boundary, even small avalanches may be effected by the 
finite system size. However, given the nature of these di- 
rected model's stationary states and the fact that all sites 
in the lattice are equivalent to the extent that the sys- 
tem is infinite, one may reduce edge effects by introducing 
sand at a fixed site a maximum distance from the bound- 
aries in randomly constructed stable states. Mass data 
from the simulations was obtained by seeding avalanches 



at the top left lattice site for models I, II, IV, VI, and 
VII and in the middle of the top row for for models HI, 
V, and VIII. This mass was histogrammed, logarithmi- 
cally binned and fitted to distributions of the form 'P(s)oc 
s^~'^ where V{s) is the probability of the occurence of an 
avalanche with s toppled sites. As can be seen from the 
plots in Fig.||, in which model V was used, the deviation 
from power law behavior in the simulation where sand 
is dropped on a fixed site is much more marked than in 
the simulation with a random drop site. This eliminates 
much of the guess work in finding a suitable endpoint 
for the linefitting routine. Notice also that the region of 
power law behavior in the run with the fixed drop site 
is at least as large as that of the run with the random 
drops even though it has a considerably smaller lattice. 




Avalanche Mass 



FIG. 2. The lower curve is calculated for model V by seed- 
ing at a fixed site and using a random background on a 
500x1000 lattice, the upper curve by seeding at random lo- 
cations on a 3000x7000 lattice. The exponent in this case is 
2.43 

For the different models we found a range of values for 
T from 2.33 to 2.50. The data is summarized in Fig||, 
and exponents are given in Fig. |l]. Although a general 
trend of increasing r with increasing was observed, it 
was also found that t varied with changes in Mi for a 
fixed he- 

In model II, avalanches occur with the same distribu- 
tion as the exactly solvable model I. This is understand- 
able by observing that neither model allows untoppled 
sites, or "holes", in the interiors of their avalanche's con- 
stant time surfaces. It was this fact which was essential 
to DR's mapping of the evolution of the perimeter of an 
avalanche onto the problem of annhilating random walk- 
ers. Model II, may be similarly mapped, although with 
a different distribution over the walker's step sizes. The 
difference in possible step sizes affects the variance of the 
walk, but not the exponent for the time of intersection. 
The distributions for the lifespans of these walkers are 
then the same as that of the avalanches and are identi- 
cal for both models. However, it is possible for another 
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choice of Ni with /ic = 3 to produce avalanches with holes 
and complicated perimeters not suitable for description 
by random walkers. Model III is one such case as can be 
seen from Fig.^, and indeed, a different value for r(see 
Fig.|l|.) is obtained from its simulation. In general, as he 
increases and the toppling vectors change, different pat- 
terns of holes are allowed, and different exponents result. 
This suggests that any classification of directed models 
must include a description of hole formation. 




Avalanche Mass 

FIG. 3. Data for the models described in Fig 1. The small- 
est exponent is from Models I and II, the largest from a sim- 
ulation of an approximation to the infinite range model de- 
scribed in the text with an N of 100 and v of 2. The data has 
been normalized to make the frequency in the first bin the 
same for all datasets. 

An analytically tractable model has been devised to 
explore the limiting case of infinite he- The model lives 
on a 2d MxN lattice with he = N + v — l where v is some 
positive integer. A site becomes unstable if it possesses 
N -\-i> OT more grains. The ith unstable site is relaxed by 
subtracting the toppling vector Ti. The components of 
Ti are given below. 



iV + 1/ if s = i 
\Ti\s = < — 1 for s G the next row 

s ^ i and s ^ the next row 



(2) 



In words, if a site becomes unstable and topples, N + v 
grains of its sand are removed and one grain is added to 
every site in the next row so that v grains are lost with 
each topple. Labeling the lattice sites by the convention 
which calls the upper left lattice site one and numbers 
sites across the first row to TV beginning the next row 
with TV -|- 1 and so on, one finds the toppling matrix to 
be lower triangular. As stated earlier this implies that the 
critical state consists of all stable configurations placing 
this model in the same category as the rest of the models 
considered in this treatment. 



FIG. 4. This is an avalanche for model III chosen to il- 
lustrate the hole formation. Note the presence of untoppled 
sites(white) within the boundaries of the avalanche 

Note that with this model and for an avalanche be- 
ginning with a single unstable site, the constant time 
surfaces and therefore toppling activity are restricted to 
a single row. This allows the row by row evolution of 
an avalanches' mass to be described by the properties of 
a stationary Markov chain with the row number serving 
also as the time parameter. 

Let Sn be the random variable associated with the 
number of sites which topple on the n*'' row. Since a sites 
height in the stationary state can take values from to 
iV -|- — 1 with equal probability, the conditional prob- 
ability that a single stable site has z grains and topples 
given that i topple on the previous row is 



V[Sn+i =1|5„ ^i]^ V[z> N + i^-i] 



(3) 



The probabilities of different sites toppling on the same 
row are independent so the transition probabilty, p^ , of 
j sites toppling given that i toppled on the previous row 
can be written as, 



VlSn+i = j\Sn = i] 



N - 



N - 



N - 



N-j 



(4) 



It is well known that a binomial distribution of this form 
approaches a Poisson distribution in the limit of large 
N so that the above transition probability may be reex- 
pressed as, 



N 



lini Pi J 



(5) 



Letting T be the random variable associated with the 
avalanche lifetime, (or equivalently the number of rows 
which have at least one toppling event), it is clear that 
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P[T <n]= P[5„+i = 0] (6) 

and so 

P[T < n] = P[Si = 0] + P[Sn+i = 0\Si = l]P[Si = 1] 

(7) 

for avalanches begun by adding a single grain to the first 
row. P[Si = 0] and P[Si = 1] are constants for a fixed 
system size and P[Sn+i = 0|iSi = f] can be expanded in 
terms of the transition probabilities as, 

N N N 

-P[5„+i = 0|5i = f] = ^ ^ . . . ^ Ps„oPs„_is„ • ■ -Pis 
s„=os„_i=o S2=a 

(8) 

. Carrying out the above sums for different n's leads to 
the following sequence of conditional probabilities. 

P[Si = 0\Si = 1] = 
P[S2 = 0\Si = 1]= e"i 
P[S3 = 0|5i = f] = e^['5==o|5i=il-i 

P[5„+i = 0|5i - f] = e^'[5.=o|5,=i]-i 

Expanding P[5„ = 0|iSi = 1] in powers of ^ and using 
the above recursion to find the cocfhcients, one finds that 
to first order, 

P[Sn = 0\Si = f] = 1 - - (fO) 

n 

. So that for long lived avalanches in large systems, we 
find that the probability of an avalanche lasting longer 
than n time steps is 

P[T>n]~— (ff) 

where A = 2P[Si — 1] and a = 1. The average flux of 
grains in and out of a row in the steady state must be 
constant. This implies that the probability of the mass 
of an avalanche being greater than m is proportional to 
m?~'^ where t is given by the relationship r = 2 + a/(l + 
a) This value of 2.5 for t matches the value obtained 
by numerical simulation to within ±.0f for N — 100 and 
v=2. 

From the point of view of the analogy with equilibrium 
critical systems, the continuous variation of the avalanche 
exponent with the range and details of the interaction 
is not understandable. There are of course equilibrium 
models with varying exponents, such as the Ashkin- Teller 
model, but we see no connection of the SOC models we 
have considered with these models. Varying exponents 
have also been obtained in dissipative SOC models, such 
as the one introduced by Olami, Feder, and Christensen 

, but in these models, the dissipation, which effectively 



links every site to the boundary, is a long ranged interac- 
tion, and so the variation of exponents is at least not in 
contradiction with the intuition gained from equilibrium 
systems. [|| 

However, in our case, the models intermediate between 
the DR model and the mean field model we have intro- 
duced are all short ranged in the usual sense of equi- 
librium critical phenomena. They should all have the 
same critical exponents, which is Ben-Hur and Biham's 
conjecture, and they do not. We do not think that the 
natural point of view from this perspective, that we are 
just not yet in the critical region for the data we have 
shown, and they will all eventually cross over to a com- 
mon exponent of 2.33, is tenable. By using the fact that 
the recurrent set is the stable cube for our models, we 
eliminate the problem of long transients or uncertainty 
in the exponent due to edge effects. For model three, 
the system size used is about 100 times the interaction 
range, which is the only scale in the problem. If there was 
going to be a cross-over, we should have seen it. We con- 
clude that there is no universality in those 2-d directed 
models whose avalanches from recurrent states contain 
holes. The extent to which the breakdown of universal- 
ity extends to other SOC models remains to be seen [ p^ . 
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